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Abstract. By instantaneously changing a global parameter in an extended 
quantum system, an initially equilibrated state will afterwards undergo a com- 
plex non-equilibrium unitary evolution whose description is extremely chal- 
lenging. A non-perturbative method giving a controlled error in the long time 
limit remained highly desirable to understand general features of the quench 
induced quantum dynamics. In this paper we show how integrability (via the 
algebraic Bethe ansatz) gives one numerical access, in a nearly exact manner, 
to the dynamics resulting from a global interaction quench of an ensemble 
of fermions with pairing interactions (Richardson's model). This possibility is 
deeply linked to the specific structure of this particular integrable model which 
gives simple expressions for the scalar product of eigenstates of two different 
Hamiltonians. We show how, despite the fact that a sudden quench can create 
excitations at any frequency, a drastic truncation of the Hilbert space can be 
carried out therefore allowing access to large systems. The small truncation 
error which results does not change with time and consequently the method 
grants access to a controlled description of the long time behavior which is a 
hard to reach limit with other numerical approaches. 



1. Statement of the problem 

Due to fast dissipation rates, the unitary evolution of out-of-equilibriuni quan- 
tum systems remained for a long time a question which remained mostly of academic 
interest. However, recent developments in the field of cold atomic systems [1 , in 
which tunability of the Hamiltonian combined to very low dissipation made this 
problem experimentally accessible and created a strong interest O [3l |4l |5] in try- 
ing to draw a general theoretical understanding of far from equilibrium quantum 
dynamics. The interaction quench is a simple way to create a prepared out-of- 
equilibrium quantum state of the system. It is realized by supposing a quantum 
system is initially in an eigenstate (the ground state) of Hamiltonian Hg^ . Here go 
is some global tunable parameter of the Hamiltonian. At time t = 0, the parameter 
go is instantaneously changed to g and as a consequence the initial state IV'go) 
longer an eigenstate of Hg. This process puts the system in a far from equilibrium 
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state, whose subsequent evolution is fully determined by the new Hamiltonian Hg 
and the initial conditions. 

The resulting dynamics of the wave function are easily formulated using elemen- 
tary quantum mechanics since they simply are given by solving the time dependent 
Schrodinger equation: 



(1.1) |^(t))=e-^.*|^M). 

Naturally, the evolution is greatly simplified by projecting the system onto the 
eigenbasis of Hg defined by a set of eigenvectors |V'g) with eigenvalues uj" . 



(1-2) \m)^T.^''^"'(^'9\^9o)\rg)- 

V 

The resulting time-evolved average of any given operator O can also be written 
as a sum over form factors weighted by the projections of the initial ground-state 
onto the eigenbasis of the final Hamiltonian: 



(1.3) m)\om)) ^Y^e'^^-"--" [W9o)ii^^\rg')\ 

Not only does one need to be able to compute eigenstates and eigenenergies 
of both Hamiltonians, the exact treatment also requires the capacity to compute 
overlaps between the eigenstates of two different Hamiltonians. We will show how 
this apparent hurdle is simply dealt with, when using the results from the algebraic 
Bethe ansatz (ABA) within the fermionic pairing model treated here. 

Moreover, as is readily seen, obtaining exact results necessitates a double sum 
over the full Hilbert space. Since its dimension generally grows factorially with 
system size, this is apparently heavily limiting as far as reachable system sizes is 
concerned. In the coming sections, we also show how the quench populates predom- 
inantly a small subset of the Hilbert space therefore allowing a drastic truncation 
with minimal loss of information ultimately granting access to systems much larger 
than one could naively expect. 

Most of the results presented here were previously published in ^ . This paper 
however supplements this reference with additional details on the truncation scheme 
and also presents novel ideas on the relationship between weak and strong coupling 
eigenstates of this particular system. 

2. System 

We treat a model of spin 1/2 fermions in a set of single particle energy levels 
which can accommodate up to two electrons (of opposite spin). Electrons interact 
through a typical Bardeen-Cooper-Schrieffer (BCS) uniform s-wave pairing term 
which retains coupling only between time-reversed states: 



a a a, 13 

The model is basically a discrete version of the celebrated BCS Hamiltonian 
to which it reduces in the thermodynamic limit An Anderson pseudo-spin 
representation [8] also exists by defining — Cq^|Cq_|, which leads to 
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(2.2) HBcs{{ej},g)^X.'^^!-9 E ^a^^' 

j = l a,/3=l 

Here the N levels are the unblocked levels, i.e. the ones which are not occupied 
by a single electron and therefore participate to the dynamics [9j. This model, 
although initially developed by Richardson [9j to describe pairing in nuclei, has 
also been used more recently to describe the physics of ultrasmall superconducting 
metallic grains '10'. 

Irrespective of the distribution of the energy levels and of their degeneracies, 
the model is integrable and can be diagonalized fully using the ABA. 

2.1. The eigenstates of the systems. Through the use of the algebraic 
Bethe ansatz |10|, \12\ I13j . one can access a very compact representation of the 



eigenstates of Hamiltonian 2.2 For any given set of M complex parameters Wj we 
shall call rapidities, the repeated action of the C{wk) operator of the monodromy 
matrix on a fully down polarized state, builds a "general" state: 

j\/ M / N c+ \ 

(2.3) \{w,}) = n C{w,)\ Hi ... i) ^ n E I - ^) • 

fe=l k=l \Q = 1 / 

All of these M rapidities states have a definite ^-projection of the total spin 
operator which is simply given by M — y. For any value of g, the unnormalized 
eigenstates of the BCS Hamiltonian are all built in this fashion. For a given cou- 
pling constant g, the complete set of eigenstates correspond to the general states 
constructed out every set of rapidities which are solution to the Richardson (Bethe) 
non-linear algebraic system of equations: 

(2.4) J^1,...,M. 

The eigenenergy of one of these states is given (up to a constant) by the sum 

M 

(2.5) Ei{wj})^Y.'^3- 

At zero coupling, the full set of fixed magnetization (fixed M) eigenstates of 

N 

Hamiltonian _ffo({£j}) = ^i^j known to be given by the various colinear spins 

states with M spins aligned in the +z direction and the remaining ones in the ~z 
direction. Each and every one of these states can be expressed in the form |2.3| by 
using a set of M real rapidities Wj . By giving the M rapidities the values of the M 
energy levels at which b, + ^ spin is found {{wj} = {e*}), this leads to 

(2.6) IK})(x l^n^+j UU...i), 

which is just the state with AI flipped spins at the chosen levels {e*}. 
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In the opposite regime of infinitely large coupling, the Hamiltonian reduces to 



N 

(2.7) H^{{e,}) - -5 E ^^^^ = -gSLS^ot = "5 [J^ " (J^^ + ^1 , 

a, (3=1 

where J is the total angular moment operator of the N spins system and is its 
z component. The eigenvectors of the Hamiltonian are therefore known through 
the procedure of angular momentum coupling of N spins | leading to sets of eigen- 
states of and J^. As mentioned previously, the z-projection of the total angular 
momentum is equal to the number of rapidities M. For a fixed value of M, one can 
absorb the Af-dependent terms in a global shift of the energies and consequently the 
eigenspectrum consists of i^f^_My.{My. states subdivided into degenerate sub-bands 
of total angular momentum j with eigenenergies, and degeneracies given by: 

(2.8) E{j) = -<?[j(j + l)] 

7V!(2j + 1) 



= (f+., + l)!(f-.): 

In the strong coupling limit, the solutions to Richardson's equations naturally 
defines an alternative representation for these eigenstates. It is known that for a 
state defined by M rapidities, the various solutions will be defined by a number 
r of divergent rapidities. These rapidities can be expressed in the large g limit as 
Wi = Lig, where the L^s are the zeros of the Laguerre polynomial 



(2.10) L;1-(™-2M)(^^)^0 y^^^ 

There is a direct correspondence between the number r of diverging rapidities 
and the total spin eigenvalue j. Looking simply at the energies, it was shown |14j 
that 



(2.11) r^j + M-- 

The remaining roots (which remain finite) defining the various eigenstates can 
then be found by solving the following set of simplified Richardson's equations 
which involve only the non-divergent M — r rapidities: 



N M — r 

(2.12) = j = l,...,M-r. 

This particular set of equations still requires a numerical approach in order to 
find every possible solution. 

One must realize that within a given degenerate subset both sets of eigenstates 
(Bethe Ansatz and total spin eigenbasis) are chosen differently. This is of course 
not an issue in the g = oo case, but the Bethe ansatz representation compares 
advantageously to the total spin eigenbasis when trying to expand in powers of 
^. Whereas in the former, one can directly expand the Richardson equations in 

powers of ^ in order to describe the system |14] . in the latter representation we 
would need to apply perturbation theory on a set of highly degenerate levels which 
would therefore involve diagonalization of large matrices. 
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2.2. Numerical solution. Except in the previously discussed extreme cases 
of zero or infinite coupling, no analytical solution is known to the Richardson equa- 
tions and one therefore needs to resort to numerical techniques in order to study this 
particular system. Many efforts have been made to construct efliicient algorithms 
for solving Richardson's equations |15L 116] but it still remains a difficult task to 
achieve stable and fast computation of the various solutions. As is always the case 
when dealing with non-linear systems of algebraic equations, the convergence of 
any algorithm towards a given solution highly depends on the capacity one has to 
generate an appropriate initial approximation of it. 

In the absence of any known general features of the solutions in the intermediate 
coupling regime, it makes it quite natural to use a scanning procedure starting from 
the known solutions at g = and increasing the coupling in small steps in order 
to reach solutions at a given value of g. This guarantees an adequate trial solution 
from which an iterative procedure can be carried out successfully. 

The solutions are such that rapidities are either real or form a complex conju- 
gate pair (CCP) with another one. At the critical values of g at which two real ra- 
pidities collapse into a CCP (or vice- versa), Richardson's equations expressed in the 
form |2.4| have cancelling divergent terms which require appropriate modifications 
in order to be dealt with numerically. These forming and splitting apart of CCPs 
happen exclusively when two rapidities, say wi and 11:2 are equal to one another 
while at the same time being equal to one of the energy levels Cc G {ci, e2---ejv}- 
This allows a rewriting of the equations using the two real variables: 



(2.13) 



A-|_ = 2ec — wi ~ W2 
A_ = (wi - ^2)^- 



Contrarily to wi and W2, these variables maintain a finite derivative with re- 
spect to g while going from real to complex. One can additionally transform the 
Richardson equations for wi and W2 in the following two equations: 



(2.14) 
(2.15) 

where 



4A++((A+)2-A_)Gi 
((A+)2~A_)A_G2 







4A+ = 



9 ^ 

N r 

(2.17) G2 = -J2 



1 



1 



Wi - ta W2- £c 
1 



{wi - ea)(w2 - £«) 



M 
M 



1 



1 



Wi-Wk W2-Wk_ 
1 

[wi - Wk){w2 - Wk) 



Rapidities which form a CCP at a given g can split apart into two real rapidities 
at a larger coupling and can also reform a CCP with a different rapidity.This means 
that one might need to redefine a new set of variables (and equations), at every step 
in g. This is simply achieved by "pairing" any two real rapidities Wi, Wj which have 
the same closest single particle energy level ec, and using for them the variables |2.13| 
For these rapidities we also use the modified version of the Richardson equations 
written above. This always results in a system of AI equations depending on M 
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variables, but it can be modified at every step in g in order to adapt to the local 
structure of the solutions. 

Independently of the pairing used at a given point, the Jacobian is always ob- 
tainable analytically in terms of these variables and it therefore allows a straight- 
forward use of Newton's algorithm for non-linear algebraic equations. Despite the 
lack of guaranteed stability (unavoidable in a general set of algebraic equations), 
one can still manage to compute rapidly any solution in this way, provided we ad- 
just the steps in g in such a way that at every g we correctly evaluate the rapidities 
which need to be paired. In a case where numerous pairings- unpairings happen on 
a small g interval, one simply needs to reduce the steps in g in order to make sure 
every rapidity crosses only one critical point per step. Since the scanning process 
gives access to a local evaluation of the derivatives of rapidities, we can always get 
a good evaluation of the desired Ag for the coming step. 

Figure [T] shows the real part of rapidities for a few solutions of the Richardson 
equations (examples showing also the imaginary part can be seen in other papers 
[l7j for example). These examples show how one can conjecture, for a given g — Q 
solution, the number of diverging rapidities that this state will have when deformed 
from to infinite coupling. 




Figure 1. A few examples of numerical solutions for M=8,N=:16 
(real part only) showing how to relate the 5 = structure to the 
one &t g — oo. We respectively get 2,3 and 5 divergent rapidities 

One can understand the passage of a given state from its uncoupled form to 
the resulting large coupling state, by looking at the contiguous blocks of occupied 
and unoccupied states at g = 0. Any contiguous block of m unoccupied states 
will refrain up to m rapidities above it from diverging to infinity. Starting from 
high energies, one looks at the first block of occupied states, and the first following 
block of empty states. Supposing they respectively contain pi states and hi states, 
if hi > pi, the pi rapidities will remain finite. If hi < pi then hi of these rapidities 
will remain finite and the remaining pi — hi can be added to the following block of 
occupied states. Supposing the next blocks of occupied (empty) states contain p2 

{P2 
P2+P1- hi 

depending on what happened in the previous block. In the three examples presented 
in Figure[2| we respectively find for the leftmost case: {pi = 1, /ii = 2); {p2 = 2,h2 = 
1); {p3 = 2, /13 = 1); {p4 = 3, /14 = 3) leading to two diverging rapidities (here they 
form a CCP), in the middle panel we have {pi — 3,hi — 1); {p2 — 2, /i2 = 1); {ps — 
l,/i3 = 1); (p4 = 2, hi = 2) giving 3 diverging rapidities (two of them forming a 
CCP) and the last case has {pi = 3, hi — 3); {p2 = 5, /12 = 0) giving 5 divergent 
rapidities (including 2 CCPs). Although giving results equivalent to the algorithm 
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proposed in |16j . this approach can also provide some insight on approximative 
values of the finite rapidities in the large g limit. This could give access, in this 
limit, to initial guesses close enough to the real solutions to allow direct computation 
of the various solutions to eq. |2.12[ 

This procedure seems to work well for any state which is half-filled or less, 
since in that case the diverging rapidities all have real parts going towards —oo, the 
Laguerre polynomials (2.10 1 involved having all their zeros in the left half plane. At 
higher filling factors this is no longer true, but the alternative choice of representing 
the eigenstates as the repeated application of the operator B{wi) on the fully up- 
polarized pseudo-vacuum: IlfcLi {j2a=i wk-e ) I ^1^1^ ■'■ ^o^^l'^ probably allow 
similar arguments to be used in this case. 



3. Quench Matrix 



Having shown how one can access individually any eigenstate of the system, we 
now need to be able to realize explicit calculations of the projection of the initial 
state of the system onto the eigenbasis of the final Hamiltonian. For a general 
Hamiltonian and even for a general integrable Hamiltonian this problem requires 
full diagonalization, in a given basis, of the initial and final Hamiltonians. The 
factorial growth of the Hilbert space dimension (and therefore of the Hamiltonian 
matrix) limits one to very small system sizes. 

For ABA integrable systems there exists a formula, due to Slavnov |18j . which 
gives an expression of the overlap between an eigenstate of the system and a state 



built identically (eq. 2.3 1 but for a general set of rapidities. This projection can 
be written as the determinant of an M by M square matrix, which in the model 
discussed here reads explicitly |13j : 



(3.1) {{w,}\{v,}) 



M 



M 



^BetG[{w,}]BetG[{v,}] ^ 

[[{Wa- Wb) [[{Va- Vb) 

a<b a>b 



Bet J, 



with the matrix elements of J given by: 



(3.2) Jab 



Vb - Wb 
Vb - Wa 



N 

E 



2E 



1 



(Wa ~ Ulc)iVa - Wc) 



k—1 c^a 

The M by M Gaudin matrices G, whose presence is needed to normalize both 



eigenstates, are themselves obtained through the limit {wj} — *■ {vj} of eq. 
which gives: 



3.1 



(3.3) Gab[{Vj}] 



y ^ 2V ^ 

2 

(Va - Vb)^ 



(Va - VcY 



a — b , 



a^b. 
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In the precise case treated here, the C{w) operators have no exphcit dependency 



on g (see eq. 2.3 1. Consequently eigenstates of the Hamiltonian at any value of g 



are all general states built a la eq. |2.3[ This allows one to use formula |3.1| directly 
to compute the relevant overlaps between two sets of eigenstates {iljg\4>gg)- One 
must understand that a different model might not share this property and if a 
determinant representation exists for these overlaps in a general integrable model, 
it remains to be found. 



3.1. Weak coupling to strong coupling. For M — 8 rapidities and = 16 
energy level, the dimension of the Hilbert space is 12870 which is sufficiently small 
to allow us to compute the full set of eigenstates. Using Slavnov's formula |3.1| we 
calculate the overlaps {ij;g\tpg^) of the go ground state with the full eigenbasis at g. 
For a few quenches, focusing on equally spaced levels (e^+i — = 1), the squared 
norm of these overlaps are plotted on figure [2] as a function of the energy of the 
state l-^p. 




Energy {in units of tiie interievel spacing) 



Figure 2. First column of the quench matrix (ground-state over- 
laps) for several quenches. In all plots N = 16, M = 8 and the 
ground state energies (represented by vertical lines) have been 
shifted for clarity. Top: Decomposition of the g ~ ground-state 
with states at g = 0.05, 0.5, 0.95. Center: Decomposition of several 
initial ground-state go = 0,0.15,0.3,0.5 in terms of the states at 
(7 = 1. Bottom: Decomposition of the go = 1 ground-state in terms 
of 5 = 0.95,0.55,0.15,0 states (from ref. [6]). 

In the limit of a small quench from g = to g = Ag, first order perturbation 
theory theory can be used. It gives non-zero overlaps with the g = Ag coupled 
ground state and the set of states which differed (at 5 = 0) from the ground state 
by a single particle-hole excitation (one up spin moved to a different level). The 
following order would allow non-zero overlaps with two particle-hole excitations 
states and so on. This fact is clearly seen in the smallest quench presented in the 
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top panel {g — 0.05) of figure [2] showing a set of decreasingly oveiiaping subbands, 
each associated with a different number of particle hole excitations 

When the quench is carried from an initial weak-coupling ground state to a 
stronger coupling (the two top panels of figure [2]) , one finds that a surprisingly 
small fraction of the final Hilbert space is importantly populated by the quench. In 
fact we see, in the distribution of energies, the large coupling band structure (see 
eq. 2.8) forming and notice that in each of the bands the projection on a single 



state of the initial ground state heavily dominates. We could easily verify that these 
states are the finite g deformations of a given set of g = configurations pictorially 
shown on figure [s] This set of states is obtained (at <? = 0) by promoting contiguous 
blocks of Np < M rapidities from right below to right above the Fermi level (FL). 
According to the procedure outlined in section |2.2[ they all lead, a,t g ^ oo, to 
states whose number of rapidities remaining finite is equal to Np. In fact, they 
are the lowest energy states having a given number of particle-hole excitations (at 
g = 0) equal to the number of finite rapidities (at g — > oo). 



• • • • o o 

• • • o o o 

• • o o o o 

• o o o o o 

o o o o o o 



o o o o o o 

• o o o o o 

• • o o o o 
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• • • • o o 

• • • • • o 



Figure 3. Pictorial representation of the "single block states" ob- 
tained by promoting contiguous blocks of A'p < M rapidities from 
right below to right above the Fermi level. 

For any quench from go = to g E [0,1], this set of A/ + 1 states always 
account for more than 60% of the total amplitude KV'gl^")!^ of the wave func- 

tion expressed in the final eigenbasis. Figure |4] shows this total contribution when 
quenching from an uncoupled initial state. Despite the presence of a dip in their 
contribution when quenching to the intermediate coupling regime, they clearly re- 
main strongly dominant for finite size systems. 

It is then natural, when looking for additional large contributions, to assume 
that states built by slightly deforming the preceding set would be dominant. We 
therefore add to the reduced eigenbasis the set of states built at 5 = by adding 
a single particle-hole excitation to the block-states. Figure [5] shows, for one of the 
block states discussed previously, the g = configurations created by adding this 
single excitation (either above of below the Fermi level) . 

One can then add a set of states built with two particle-holes excitations over 
the block states. We include every state obtained with one particle excitation above 
the Fermi level and one hole excitation below it. These would be constructed by the 
tensor product of the possible states above the FL, as seen in the first set of states 
in Figure [5j with the possible states below the FL in the second set. Naturally, 
states with two excitations above (or below) the FL are also possible. A fairly 
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Figure 4. Total contribution of these states to the amphtude of 
the initial state at go = 0, as a, function of interaction (from ref. 

m)- 
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Figure 5. Pictorial representation of the single-excitations for a 
given block state. One should understand that the same process is 
applied to every one of the block states (see Fig. [s]) 



small set of these is actually included, since the study of the quench matrix made 
for a smaller system indicates that they should be dominant. The treated cases 
are built by promoting any two rapidities (absence of rapidities) above (below) the 
Fermi level to any contiguous set of two free (occupied) level above the FL as shown 
in figure |6] 

In the case of M = 16, = 32, the inclusion of all the aforementioned states 
proves sufficient to obtain more than 97% of the weight of the wavefunction after 
the quench, i.e.: l('0g|V'°)P > 0.97. This figure holds for any final value of g 

V 

between and 1. One must understand that for a given value of the final coupling 
only a part of these states will have an important contribution to the wavefunction. 
As one can readily see in Figure |2] the weight is shifted to different parts of the 
Hilbert space as g changes. For example, when doing a large quench the overlap of 
the final ground state with the initial one will become quite small, but for a small 
quench it would be the dominantly populated state. It is therefore possible to get 
an accurate projection of the uncoupled ground-state onto the coupled eigenbasis 
by keeping less than a 1000 states (out of the 601080390-dimensional Hilbert space 
for M — 16, — 32) for a given quench. However, in order to be able to treat 
the spectrum of quenches discussed here we need to compute only a total of 7675 
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Figure 6. Pictorial representation of the doubly-excited states 
which are kept for a given block state. The same process is applied 
to every one of the block states (see Fig. |3]). Moreover, one should 
remember the inclusion of states featuring one excitation above 
and one below the FL (not represented here but discussed in the 
text). 



states. This drastic reduction of the effective Hilbert space is the key element which 
allows us to treat relatively larger system sizes with a minimal error. 

Since this truncation scheme is carried out in the eigenbasis of the Hamiltonian 
driving the time evolution of the system, the error will not change with time. In fact, 
the time evolution only affects the phase of various components of the wave function 
and not their amplitude. Apart from the initial condition, the wavefunction (and 
by extension expectation value of operators) at any time does not depend on the 
previous history. For very large times, it can therefore be calculated directly and 
with a precision set by the small truncation error. 



3.2. Strong coupling to weak coupling. Quenching form a strongly cou- 
pled initial ground state down to a weakly coupled system brings a radically different 
structure. Looking at figure |2] we see that although a small quench still brings a 
quench matrix naturally strongly peaked at low energies (with a roughly exponen- 
tial decay in energy), when the difference between initial and final g is made larger, 
the distribution flattens significantly. 

One can understand this tendency by noting that the infinite coupling ground 



state is characterized by having every rapidity present diverging (see eq. 2.5 1. In 
this case we have hm^j^oo C{v) oc ^ ■ Sf . The repeated action of the operator for a 
set of diverging rapidities therefore make the M rapidities ground state a uniform 
superposition of every g = A/-flipped spins state, i.e. 



M 



(3.4) K^oo>«(^E^g UU..4). 

The g oo ground state therefore has a uniform projection onto the g ^ eigenba- 
sis. In this extreme scenario, any truncation of the final Hilbert space is impossible. 
However, when starting from a finite albeit large initial coupling, it appears that 
the nearly monotonic energy behaviour could allow one to reduce the effective di- 
mension of the problem by removing the largest energy states. The slow decrease 
of the overlaps in energy would still make the problem of very large to very weak g 
quench hard to treat using this technique. 
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4. Form factors 

Being in a position to compute the eigenstates, eigenvalues and the quench 
matrix elements is sufficient to access the time-evolved wave function of the system. 
However, to compute the average value of an observable O, one still needs tractable 
expressions for the form factors: 



(4.1) mo\i^^}. 

of the operator we are interested in. In this work we will look at time evolution of 
the off-diagonal order parameter defined as 



a/3 



In equilibrium and in the thermodynamic limit, this quantity is closely related 
to the BCS gap making it a suitable parameter to quantify the superconducting 
tendency for finite systems even out of equilibrium. One could also use in a similar 

fashion the canonical order parameter J2a \J \ ~' i^a)^ (as in |19L I20jl. but since 
both quantities showed the same qualitative behavior we only present the former 
here. 

Although in the infinite couphng case we have ^(^g 5^ IV'g ) ^ ^ti,i^7 since 

a/3 

^ ^ Sg CX Hryo, at any finite coupling it requires computation, for every couple 

a/3 

a, of the form factors {1p''g\S+Sg\i;^^). 

Using Slavnov's formula, the solution to the inverse problem, i.e: 



S+ = lim(u-ei)C(u) 

(4.3) = YiTn{u-€,)B{u), 

u — >ei 

and determinant properties, it was shown in f6] that these form factors for S~Sp 
(and consequently S~S'^ too) can be written as a sum of M determinants. This 
was done by generalizing to the case {v} ^ {w} the method of Ref. |1_7J, starting 
from the double sums given in Ref. [13j. For a ^ P we have 

n(^)n(-^-.) M 

(4.4) (M|5-^;iW) = ^ / -E?^det:^.%' 
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where, defining Aab = JabY\_i'"c ^ ^b) {Jab being defined by eq. 

c^b 

elements are given by 



3.2 1, the matrix 



(4.5) 



O aq 

jq 



I\k^b,qi'^'' - Wb) Wb-ea 



Uk^b+l,qi'>"k - Wb+l) Wb+1 - 



A, 



6+1 1 



b<q-l, 



Aa„-l+2 



{Wq - ei3)(Wq-l - Eq) 



n 



Y\ {Wk — Wq-l) 



Wq-l — Wq 



Vc - (.(} 
Wc — 6/3 



k^q-1 
l/(Ua — ^af, 

Aab, b>q. 



{Va - ta)'^iVa - tpY ' 



For a = (3 \t can be easily related to the form factors for S*^ which were already 
published elsewhere |13j . 



5. Order parameter evolution 

Having a numerically tractable expression for every necessary element, is it 
straightforward to compute the time evolution of the off-diagonal order parameter 
(eq. 4.2 1. The upper left corner of figure It] compares the time-averaged value i.e.: 



(5.1) * = ^lini^ r dt^{t) = l«lOl'Kli?7E^^^^"I^P 

of the off diagonal order parameter with the thermodynamic limit given by the 
mean-field BCS theory. The non-equilibrium 'phase diagram' for iV — > oo shows 
that the final asymptotic value of the canonical gap Aco defines a universal curve 
when expressing Aoo/Ag vs Ag„/Ag, where Ag is the equihbrium value (ground 
state expectation value) at coupling g [19] . In term of the off-diagonal order 

parameter, the gaps are given by A^ = .9^((V'gl Ea,/3=i -S'^-S'^IV'g)/^^ " 1)/M HU 

and we use A^c, = g\J /M — l) /M. One can clearly see that for Ag^/Ag > 1, 
which corresponds to a quench going from large coupling to weak coupling, the 
results show stronger deviations from the BCS mean-field result obtained in the 
limit N oo. This is to be expected, since for a finite system, at weak coupling, the 
superconducting correlations between electron pairs are suppressed when compared 
to the infinite system BCS result. In the limit g 0+, the infinite system still shows 
an instability which leads to a BCS wavefunction, but for any finite system size, 
the small coupling limit leads to a ground state formed by uncorrelated electron 
pairs and the mean-field BCS treatment loses validity. Nevertheless, in the regime 
of relatively small quenches (A^^/Ag ~ 1) one still sees that the impact of quantum 
fluctuations (which are not present in the BCS mean-field tratment) on the long 
time average gap is strongly suppressed even for the relatively small system sizes 
treated here. 

The two additional panels in figure |7] show the time-evolution (and its Fourier 
transform) of the off-diagonal order parameter. Those results are given for M — 
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Figure 7. Bottom left: OfF-diagonal order parameter evolution 
for iV = 32, M — 16. Right: Fourier transform, the various 
plots are shifted on the vertical axis for clarity. Top Left: Non- 
equilibrium finite-size "phase diagram" resulting from the time- 
averaged canonical gap obtained from the off-diagonal order pa- 
rameter, as explained in the text (from ref. [6]). 



16, = 32 which therefore necessitates the use of the previously described trun- 
cated eigenbasis. We focus on quenches from the g ^ uncorrelated ground-state 
to various finite values of the coupling. 

For general quenches, the mean- field treatment, is known to lead to integrable 
classical dynamics, for which, in the steady-state reached for quenches from small 
to large coupling, the time evolution of the BCS gap is given by the Jacobi elliptic 
function [19, 2Q\ : 



(5.2) A{t) = A+dn[A+{t-To),k], fc = 1 - A^/A 



with parameters A+ Ag and A_ for a weakly coupled g O^initial state. 
This gives rise to non-harmonic persistent oscillations which are periodic in time 
with a period given by the complete elliptic integral of the first kind: 



(5.3) T: 



^+ Jo ^1 - fc2 sin^ (f) 

The mean-field solution should therefore show a Fourier transform made of 
equally spaced peaks which, as one can see, is quite different from the results ob- 
tained here. For small values of the final g, this difference is no surprise since the 
mean-field treatment assumes a BCS-like wave function, which, for a finite size 
system in sufficiently weak coupling is not realized. However, when g ^ g* — 
(21n7V)^^, it was shown |17l |21, that the equilibrium static correlation functions 
are undistinguishable from the BCS correlations. For the largest final g values 
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shown here this inequahty is respected but the quench dynamics are still qualita- 
tively different from the expected BCS result. 

One can understand this discrepancy by simply looking at the overlaps plot- 
ted on figure [2] It was clear from the start that an instantaneous quench can, in 
principle, create excitations at arbitrarily large frequencies. What our calculations 
show is that the excitation is done in a highly non-thermal way, by creating pref- 
erentially excitations in the mid-range energy spectrum. On the other hand, static 
correlations functions depend mostly on low-energy properties. Since the equis- 
paced low energy bands which characterize the mean-field BCS regime are already 
well formed a.t g > g* — (21nA^)~^ (see the energy distribution on figure |2]), it is 
therefore expected that low energy properties are quite similar to the BCS ones. 
The quench dynamics, however, by probing higher energy properties which still 
differ strongly from the BCS behaviour, will retain a visible impact of the quantum 
fluctuations that the mean-field treatment cannot capture. As it was proposed in 
other models [4], the experimental accessibility to quantum quenches could pro- 
vide a powerful spectroscopic tool to study quantum fluctuations in the pairing 
properties of fermions. 

6. Discussion 

In this work we have reported a new approach to the study the dynamics of a 
quantum system driven out of equilibrium by a global interaction quench. It uses 
integrability which provides numerically tractable ways to compute eigenstates, 
eigenenergies and form factors of various operators. We were able to exploit the 
peculiar property of the Richardson model in order to find simple expressions giving 
the overlaps of the eigenstates of the pre and post quench Hamiltonians. 

The obtained results show that the particular structure of these overlaps, also 
allows one to design a simple but very efficient truncated Hilbert space giving access 
to nearly exact results in systems too large to allow the full Hilbert space to be 
treated. 

The generalization to other integrable models such as quantum spin chains 
[22] or atomic Bose gases |23| . would be a highly desirable feat. However, in other 
general integrable models, the current lack of known tractable expressions for the 
elements of the quench matrix hinders such progress for now. Such a representation 
would allow exact calculation of quench behavior for a large variety of potential 
experiments. 
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